Management of matter waves in optical lattices by means of Feshbach resonance 
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Mean-field dynamics of Bose-Einstein condensates loaded in an optical lattice, confined by a 
parabolic potential, and subjected to change of a scattering length by means of the Feshbach reso- 
nance is considered. The system is described by the Gross-Pitaevskii (GP) equation with varying 
nonlinearity, which in a number of cases is reduced to one-dimensional (ID) perturbed nonlinear 
Schrodinger (NLS) equations, particular form of which depends on relation among parameters of 
the problem. We analytically describe adiabatic dynamics of periodic solutions of the respective 
NLS equations, provide numerical study of ID models confined by a parabolic trap, and carry out 
numerical simulations of the matter wave dynamics within the framework of the radially symmetric 
3D GP equation. Special attention is paid to processes of generation of trains of bright and dark 
matter solitons from initially periodic waves. The results of the ID approximation are compared 
with direct numerical simulation of the original 3D GP eqation, showing remarkable coincidence for 
definite regions of parameters. 

PACS numbers: 03.75.Lm, 03.75.Kk, 05.45.Yv 



I. INTRODUCTION 



During the last decade Bose-Einstein condensates 
(BECs) attracted a great deal of attention of thephys- 
ical community (see e.g. review papers 0, 01 U and 
references therein). Creation and manipulation of spa- 
tially localized and periodic matter waves have become 
the issues of particular importance, optical lattices (OLs) 
being viewed as one of the most promising tools for man- 
agement H, 1^ ^ 1^. On the other hand, it has been 
predicted theoretically and confirmed experimentally 
P, i7|| that another tool - the Feshbach resonance (FR) 
- can be effectively employed for manipulating BECs. It 
also has been suggested that by changing external mag- 
netic field (i.e. by controlling the magnitude and even 
sign of the s-wave scattering length) one can generate 
trains of solitons starting with a periodic matter wave 
create stable localized structures such as breathers, 
2-soliton states or nearly static states with nested dark 
solitons 0, compress solitary waves up to high densities 
[Toj , generate shock waves in BECs and create single 
bright and dark solitons starting with the linear oscil- 
lator eigenstates Subsequently, a natural question 
about combined effect of FR and periodic external poten- 
tial can be posed. This issue has already been addressed 
in Ref . [l3l| . where effect of variations of the scattering 
length on dynamics of a soliton of a discrete nonlinear 
Schrodinger (NLS) equation, which simulates an array of 
BECs in a deep OL within the framework of the tight- 
binding approximation |l4j . was considered. However, 
as it was shown in |l5j |. the latter approximation does 
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not reproduce all features of dynamics of the underlining 
continuous system, what requires extending considera- 
tion beyond the tight-binding limit. 



The aim of the present paper is to consider combined 
effect of FR and of an OL on dynamics of a BEC within 
the framework of reduced one-dimensional (ID) and gov- 
erned by truly 3D mean-filed models_^ More specifically, 
we extend the ideas reported in Ref. 8] to cases of BECs 
embedded in OLs and consider the effect of FR on dy- 
namics of nonlinear Bloch states resulting in creation of 
trains of bright and dark solitons. 



The problem allows rather complete analytical descrip- 
tion in a ID case. That is why we start by listing in 
Sec.mthe main ID reductions of the 3D Gross-Pitaevskii 
(GP) equation for different relations among characteris- 
tic scales of the problem. Next, in Sec. lIIII we recall some 
relevant facts of the adiabatic theory of the NLS equation 
providing some details to the earlier results of Ref. 
This will lead us to simple picture of the effect of varying 
nonlinearity on ID periodic solutions when the lattice 
constant is small enough. In Sec. II VI adiabatic approxi- 
mation is applied to describe dynamics of periodic mater 
waves loaded in a periodic potential and affected by FR. 
In Sec. we provide numerical study of ID dynamics 
of solutions describing BEC in an OL where boundary 
effects become important due to presence of a parabolic 
trap. Sec. I VII is devoted to direct numerical simulations 
of the 3D GP equation and comparison of the results with 
predictions based on ID models. The results are summa- 
rized in Conclusion. For the sake of convenience some 
relevant information on the Inverse Scattering technique 
(1ST) and on elliptic integrals are given in Appendixes. 
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II. EFFECTIVE ID EQUATIONS 

As it is customary we start with the GP equation for 
the macroscopic wave function ^'(r, t): 



2m 



A* + y(r)* + .go|*P* 



(1) 



where we use the standard notations: go — 4:Trh?as/m, 
is the s-wave scattering length, and m is the mass 
of an atom. The external potential is given by V{r) — 
Vtrapir) + Viattix), whcrc Vtrapir) = ^ (t^i'^i +cj§a;2) 
is a trap potential, rj^ = {y,z), u!± and loq are the 
transverse and axial harmonic oscillator frequencies, and 
Viatt{x) = Viatt{x + A) is a lattice potential with the lat- 
tice constant A. The wave function is normalized to the 
total number of atoms, Af: J j^'pdr = Af. 

We consider a cigar-shaped t rap whe re the transverse 
linear oscillator length, a±_ — y^ll/rnIJ±, is much smaller 
than the longitudinal one, Oq = yjh/moj^): a±^ <^ oq. 
Then the condensate becomes a quasi-lD system. In lit- 
erature there exist several approximations reducing the 
original GP equation to different effective ID models. 
The method allowing one to do that in a controlled 
manner, i.e. appreciating the neglected terms and us- 
ing a unique well defined small parameter, is the multi- 
ple scale expansion. The resulting ID models depend- 
ing on relations among the parameters are summarized 
in the Table HI (see also j^il). To discuss them we intro- 
duce the mean healing length ^ = (87rri|as|)~^/^, where 
n ^ M / {a\aQ) is a mean particle density, and collect 
four characteristic spatial scales: {oq, a_L, A, ^}. Next we 
introduce a small parameter of the problem, e, defining 
it through the ratio between the energy of two-body in- 
teractions, ^ h^asTi/m, and the kinetic energy of the 
transverse excitations, ffi / {ma?^): 



Na 



e 



a-o 



(2) 



Then, one can single out four different cases Q. 

Let us start with the Case 1 from Table 13 for which 
the lattice period is of order of the transverse oscillator 
leng th and much less than the healing length (see also 
[l^). Introducing dimensionless independent variables 
t = ui±t/2, X = x/a^, fx = rj_/a_L, and a renormalized 
wave function ^ = ax|as|^/^\[', we rewrite Eq. in the 
dimensionless form 



Ci_ - -Ax 



£0 + 87ra|*|2U', 



Co 



(3) 



Hereafter a = sign(as), v = ujfj/uj_\_ is the aspect ratio 
of the trap, U{kx) = Viattix) /Eji is a dimensionless pe- 
riodic potential which is varying on a unit scale, = 

h?TT'^ /{2h?m) is the recoil energy, and — -£f- ^ 



Next we consider linear eigenvalue problems 

'Co'/>ng(£) = £nq4inq{S:), 



(4) 
(5) 



Hereafter n and q stand for a band number and a wave 
vector inside the first Brillouin zone (BZ), when it is rel- 
evant (cases 1 and 2 in Table P) and the index q must 
be dropped otherwise (cases 3 and 4 in Table P) . In- 
dexes li^2 stand for two transverse quantum numbers. 
When V ^ t, ipnqix) can be approximated by a Bloch 
function of the respective periodic potential. We re- 
strict the consideration to the Bloch states bordering 
the edge of the BZ where q ~ qo — tt/A, and respec- 
tively use the abbreviated notations 4>n{x) = 4>nqo{x) and 
£n = £nqo- Wc also considcr only the evolution of the 
background state {li = ^2 = 0) of the transverse compo- 
nent: ^oo(fj-) = exp(-f^/2)/0r, and £00 = 2. 

Now we look for a solution of Eq.© in a form — 
etfji + 6^-02 + • ■ ■ , where (j 



1,2, 



are functions 



of r± and of scaled independent variables ta = e"t and 
Then the leading order of the solution can be 
searched in the form of the modulated ground state 



V'l = —i=A{xi,t2)e~ 



.(£„+2).0g-fl/2^^(^^)^ (6) 



where A{xi,t2) is a slowly varying envelope amplitude 
in which by convention we indicate only the most rapid 
variables. 

The next steps are standard for the multiple-scale ex- 
pansion: substituting in © one resolves the equations 
appearing in the three lowest orders of e (see e.g. 0, 0| 
for details) and subsequently eliminates secular terms by 
imposing requirements on the amplitude A{xi,t2). For 
the case at hand this results in equation 



«5tV = -^a,V + 2a|VP^, 



(7) 



where ^ ~ \/x[2A with x = 2 |0„(x)|'*(ia;, and for 
the sake of brevity we made substitutions t2 t and 
xi X. The effective mass is obtained in the formal 
limit ao/A ^ 00: A/-^ = {<P^nq/ dq^) ^^^^ 16]. 

When the condensate size in the axial direction is not 
large enough and "boundary" effects cannot be neglected 
(Gase 2 from Tabled, the linear spectral problem must 
be modified. Since now a±/ao ^ e and thus v^x^ = 
e^iy'^Xi with 1/ — i^v and v — 0(1), the problem (PJ is to 
be redefined as follows 



/:x = -Ax+ri, Co = -dl + k^U{kx). 



Then the resulting ID equation reads 
1 



2M 



+ I' -I- 2cr|7/;| V 



(8) 



(9) 



where for the sake of simplicity iy is substituted by and 
as before {tp, t, x) stand for {y^x/2A, ^2, 2:1). 

If the period of the potential is large compared with 
the healing length (Gase 3 from Tabled, one has k = ek, 
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TABLE I: Scaling, backgrounds and effective f D equations describing matter waves in OLs. All variables in evolution equations 
are dimensionless (for the notations see the text). 



Case 


Scaling 




Background ^^(xo) 


ID NLS equation 


1 






~ Bloch function 


iVt = -^i'xx + 2cr|V'| V + i^^x^ip 


2 


ax ~ A 'N. 


~ eao 


= Bloch function 


3 


ax ~ eA ' 


- < e'ao 


^ ^-l/4g-..-5/2 




4 


ax eA ' 


■-^ e^ ~ eao 


= 1 


ixpt = -'ipxx + 2a|-!/)p^V + U{x)4; + v^x'^ip 



where k is of order one. Then k?U{kx) = e^U{xi) where 
U{xi) = k'^U{kxi) and system Q acquires the form 



(10) 



The eigenvalue problem for Co is now solved explicitly: 
(/)o(a:o) — vr"^/'* exp(— i/a;Q/2), and the final (i.e. already 
in dimensionless units) ID equation reads 



(11) 



Here U{xi) is substituted by U{x). 

Finally, if one considers not too long condensate and 
a lattice having a period of order of the mean healing 
length (Case 4 in Table P), the corresponding problem is 
reformulated now as follows 



£1 = -A, 



f2 



Lr, — —dz.. 



(12) 



The ground state as a solution of the linear eigenvalue 
problems Q becomes = 1 and the ID model reads 

idt^ = -dl^ + i^^x^^p + U{x)^j + 2cr|V'PV' (13) 

where as before the tildes are suppressed. 



III. NLS EQUATION WITH TIME DEPENDENT 
NONLINEARITY. ADIABATIC THEORY 



Cases 1 and 3 from Table allow rather complete de- 
scription within the framework of the perturbation the- 
ory, while the Cases 2 and 4 where the parabolic poten- 
tial is explicitly included in the evolution equation for the 
slowly varying amplitude, have to be studied by means 
of numerical simulations. In_the present section we re- 
produce the main results of 8] providing detail analysis 
of the asymptotics of the solutions. 



A. Stationary periodic waves 

Let us recall known formulas for periodic solutions, 
^{x,t) = %p{x + L,t), of the NLS equation (0, where 
without restriction of generality we put \M\ = 1/2 and 
define a = o-sign(Af): 



-V'rr:!; + 2ct|V|V- 



(14) 



(In the physical units the period is recovered to be a±L/e, 
and thus, in accordance with the scaling of the Case 1, 
a±L/e ~ ^ ^ A.) The simplest stationary solutions of 
such type can be obtained in a form of a standing wave 



ip{x,t) 



'u{x) 



(15) 



where u{x) is a real periodic function, u(x) — u{x + L), 
and a; is a real frequency. In this case (|14() is reduced to 
the ordinary differential equation —Uxx + + '2au^ = 
which yields 



x — Xq 



du 



(16) 



Here xq and C are the integration constants. As far as 
the stationary solution u{x) is found, a periodic wave 
moving with a constant velocity V is obtained as 



■0(2;, t) = exp 



V 



V 

, t + I — X 
4 2 



u{x-Vt). (17) 



Next steps depend on a particular choice of the pa- 
rameters C and uj. It turns out, however, that for our 
purposes another parameterization, provided by the in- 
verse scattering technique (1ST), is more convenient. 

Such parameterization is introduced in Appendix A 
and is summarized in Table ^] where the well known an- 
alytical expressions for periodic solutions are given (we 
use the standard notations for the Jacobi elliptic func- 
tions pq(a;,m) with p=c, s, d and q=n, [T7|). 



B. BEC with a varying scattering length 

As it is explained in Appendix A, controlled change of 
a shape of a periodic wave can be described by motion 
of the roots Xj on the complex plane. In particular, to 
transform a periodic wave into a soliton one has to find 
a way to move accordingly the branching points Xj. It 
turns out that the limiting transitions shown in Tabled 
and in Fig llll cannot be reached in practice because they 
require change of the boundary conditions of the problem 
and thus change of the number of atoms involved. That 
is why instead of looking for a possibility of transforming 
a periodic wave in a single soliton one can to explore a 
possibility of its transformation in a sequence (we call it 
train) of solitons equally spaced and thus preserving the 
period of the wave. This can be done by perturbing the 
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TABLE II: The simplest static periodic solutions. 



a 


periodic solution 


m 


limit m ^ 




limit m ^ 1 




1 




(6 - 6)' 
(6+6)2 


e — sm(^2;) 




-^«'*/2^^^j^l^i^j 


-1 










jye*'' *secli(r;x) 




(»7i + J72)2 




-1 


V'cn = e''* ^*^2^ cn ^ v'r/^ + ^'^x, rnj 


6 + '72 


e ' T] cosyrix) 




776*'' *secli(r;a::) 



evolution equation. This has to be done adiabatically 
in order to do not destroy the regular structure of the 
periodic wave. In particular, one can consider change of 
the effective mass by means of variation of either the laser 
field intensity or geometry of laser beams, or change of 
the nonlinearity. The latter can be achieved by means 
of the FR affecting the scattering length and controlled 
by an external magnetic filed B{t) varying with time. In 
that case a model simulating the process can be written 
as follows [la: 



as{t) = a,(0)g(i), g{t) - 



^B{0){^B{t) + Ares) 
AB{t){AB{Q) + Ares)' 



(18) 

Here as(0) is the initial value of the scattering length, 
5(0) = Bo + AB{0) is the initial field, AB{t) is a devi- 
ation from the resonant value Bq, and Ares is the width 
of the resonance. 

Assuming that dependence of g{t) on time is slow 
enough (i.e. is governed by e^t) one can repeat the argu- 
ments of Sec. m substituting go by go = go ■ g{t), what 
results in modification of each of the evolution equations 
O [respectively dl)], lEJl, dl, and In particular 

(|14|l now reads 



iilJt + i^xx~2agit)\ij\^ij^0. 



(19) 



In what follows we consider only the case when g{t) is 
a positive function, i.e. where the FR does not change 
a type of the two-body interaction. Then, in order to 
study evolution of the condensate within the framework 
of the model (|19|) we notice that substitution 



■4'ix,t) 



v{x, t) 



- p-C(t)/2 



w(x,C) 



(20) 



where 



C = C(0 = 2 / 7(i')di' = ln5(i), C(0) = (21) 

JO 



transforms (|19|l into a dissipative NLS equation 
where 

1 dg{t) 1 dC 



lit) 



2g{t) dt 2 dt' 



(22) 
(23) 



For slowly varying g{t) the right hand side of Eq. (|22(l can 
be considered as a small perturbation: \j{t)\ <C 1. It is 
to be emphasized here that slow dependence means slow 
in terms of slow time e^t, what is satisfied, for example 
by the function g depending on e^/^t. If the nonlinear- 
ity increases (jit) > 0, the case we will deal with), the 
perturbation describes growth, otherwise, when the non- 
linearity is decaying {'jit) < 0) the perturbation describes 
dissipation. 

Since g{0) = 1 it follows from Eq. (1201, that ^(a;, 0) = 
i;{x,0). 



C. Adiabatic approximation 

Under influence of the dissipative perturbation a wave 
shape is changed what can be described by variations of 
the parameters Xj assuming them to be slow functions 
of time t (the so-called adiabatic approximation). Equa- 
tions which govern their evolution can be derived by the 
following simple method 0, E] ■ 

First of all we recall that in all the examples considered 
in this section the waves are two-parametric and hence 
we have to obtain two equations of the adiabatic approx- 
imation. Since the coefficients in the dissipative model 
(|22|l arc supposed to be independent on x, the perturba- 
tion does not result in change of the period L = L(Ai, A2) 
of the nonlinear wave what can be expressed as 



d( 



L(Ai,A2) =0. 



(24) 



Here we temporarily use Ai^2 for two parameters of a 
periodic solution: thus Ai,2 = Ci,2 for the sn-wave, Ai,2 — 
r]i 2 for the dn-wave, and Ai = 77 and A2 = ^ for the cn- 
wave (see Appendix A). 

Next, it is a straightforward algebra to show that the 
integral 



N = 



'■dx 



(25) 



of a periodic solution of H22(l (below it will be referred 
to as a number of particles, assuming that there is no 
confusion with a real number of particles given by A/", see 
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Sec. evolves according to 
dC 



iV(Ai,A2). 



(26) 



Then, if the expressions for L and N in terms of Ai^2 
are known, Eqs. (|24|l and (|26l) reduce to a system of two 
first order differential equations for Ai_2- 



D. sn-wave in a BEC with OsM > 

We start with the adiabatic dynamics of a periodic 
matter wave in a BEC with a positive scattering length 
and look for a solution of (|22|l in a form of the sn-wave 
given in the first row of Table ^ Then 



L = 



8K(to) 



Af-2(a+6)[K(m)-E(m)] (27) 



where K(m) and E(m) are standard notation for the com- 
plete elliptic i nteg rals of the first and second kinds, re- 
spectively 0, I22I and m is defined through in Ta- 
ble ^1 Now, considering as functions of time, or more 
precisely of C = C(i), and substituting 1)27(1 into ((24(1 and 
((26(1 we obtain 



d^i ^ ei(ei+6)[K(m)-E(m)] 
dC 2eiK(m)-(a+6)E(m)' 
d^ ^ 6(ei+6)[K(m)-E(m)] 
dC 26K(m)-(a+6)E(m)- 



(28) 



For the sake of definiteness we assume that the pa- 
rameters a are chosen such that initially ^01 > C02 as in 
Figlllb (hereafter ^qj = ^j{t = 0)). Then one immedi- 
ately concludes that d^i/d( > 0, i.e. ^ is a monotoni- 
cally growing function of time. One can also show (see 
Appendix that ^ is a monotonically decreasing func- 
tion and thus to — )■ 1 as time goes to infinity. To describe 
the respective limiting behavior we employ asymptotics 
of the elliptic integrals I(B1(I , and reduce ((28(1 in the lead- 
ing order (when ^1 — > 00 and ^2 ^ 0) to 



d^ 

dC 



2 ■ 



d^i 
dC 



6, 



(29) 



The first of these equations is trivially solved. To solve 
the second one we observe that from ((29(1 it follows that 



da 

rfa 6 6 

This equation is also readily solved giving 



(30) 



C2 = exp(l + lna + Coa) 



Here Cq is an integration constant, which can be found 
from the (asymptotic) condition for conservation of the 
period computed from 1(27(1 : L ^ ^ In |^ ~ — 4Co. Fi- 
nally we arrive at the asymptotics: 



£.2 ~ Cie 



-Lii/4 



(31) 



Thus, one observes a scenario different from one shown 
in Figlllb.b (see the example in Fig^, below). Now 
a goes to the infinity and we obtain a "train of dark 
solitons" . This is related to the fact that the perturbation 
does not affect the period of the nonlinear wave, while in 
order to obtain a single dark soliton the period should 
go to the infinity in the process of the wave deformation. 
Meantime, the shape of the train of dark solitons indeed 
approximates to a shape of a set of static dark solitons, 
what occurs due to extremely rapid convergence of ^ to 
zero, what is illustrated in inset in Fig^. 



E. dn-wave in a BEC with asM < 

Let us now turn to the case a — —1, and consider 
deformations of the periodic dn-wave, given by the second 
row in Tabled resulting in a train of bright solitons (the 
counterpart of the transformation shown in Figs lllb .dV 
The period L and the number of particles N are now 
given by 



L 



4K(to) 



N = (?7i +ry2)E(TO). 



(32) 



Assuming that parameters rjj depend on time and sub- 
stituting (|23) in lEH) and ^ we obtain a system of the 
equations of the adiabatic approximation for 771 and 772 : 



dm 

dC 
dr]2 
dC 



i]i{m + f?2)E(TO) 



iVi + ?72)E(to) -t- (ryi - ?72)K(to) ' 

??2(??i + ?72)E(to) 

iVi + ?72)E(to) - (ryi - ?72)K(to) 



(33) 



where to is given in Tabled ( is defined by Eq. 1(21(1 . and 
for the sake of definiteness it is supposed that 772 > 771 . 

Like in the previous subsection one can show that the 
both rjj are monotonically increasing functions (see Ap- 
pendix 0. Moreover, it is evident that 772 growth more 
rapidly than 7/1. To consider the limit 1 771 — 7/2 1 0, 
i.e. TO ^ 1, we introduce rj = (771 -|- 772)72 and A77 = 
(772 — 771)72. Then using ((BID we obtain from 1(33(1 : 



dr] 

dC 



dA-q 



-A7/ln 



(34) 



dC ' A77 

Comparing these equations with 1(29(1 one observes the 
universality of the limiting behavior. Thus using l(31() 
one immediately writes down the asymptotic formulas 
(as in the previous case the constant is found from the 
conservation of the wave period) 

=^ A7y - 776-^"/'* 



(35) 



A peculiarity of the behavior of the dn-wave parame- 
ters in the limit to — *■ 1 is that both 771 and 772 tend to 
infinity (see the inset in Fig|21i, below), what is consis- 
tent with the fact that keeping the wave period to be a 
constant one cannot obtain a single bright soliton. Mean- 
time, as we could expect the difference I771 — 7/2 1 rapidly 
goes to zero, what corresponds to appearance of a "train 
of bright solitons" (see the example in Fig|3^). 
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F. cn-wave in a BEC with asM < 

Finally, we consider the case of the cn-wave (the third 
line of Table in a BEC with a negative scattering 
length. The period and the number of particles are now 
given by 

L - ^|M=, N = A^e+^^im) - (36) 

V r + T 

Substitution of these expressions into Eqs. H24(l and 126|l 

yields 

d-q _ [(772 + ^2)E(m) - ^2j^(TO)]E(m)77 

dC " ry2E2(TO) +^2[k(to) -E(to)]2 ' 

d£, _ K^K(m) - (7^2 + ^2)E(m)][K(m) - E(m)]C 

dC ~ 772E2(m) + C2[K(m) - E(m)]2 ' 

(37) 

Now "q and ^ are respectively increasing and decreasing 
functions. As in the two previous examples analytical 
treatment can be given for the limiting case, which is now 
^ ^ and ?7 — > 00. Using IjBlj) it is a straightforward 
algebra to ensure that in the case at hand we arrive at 
already familiar equations 



dr] d£, 77 

dC dC $ 



(38) 



describing the universal law of the limiting transition [c.f. 
(ED) and 123)] 



-Lri/A 



(39) 



Dynamics of respective parameters 77 and f on time is 
shown in the inset in Fig|2l3. 



IV. BEC SUBJECT TO PR IN A PERIODIC 
POTENTIAL HAVING LARGE PERIOD 

A. Stationary solution and adiabatic 
approximation 

The Case 1 from Table is the simplest one since it 
deals with the integrable NLS equation (Q. It turns 
out that the Case 3 from Table with varying nonlin- 
earity can also be treated within the framework of the 
adiabatic perturbation theory, due to the fact that some 
solutions of the NLS equation with a constant nonlinear- 
ity are known explicitly. They were considered in details 
in Ref. |23| and are partially summarized in the present 
subsection. 

To this end we notice that ansatz IjlSfl now results in 
the stationary equation 



of the NLS equation H14() . is given a priori. Indeed, one 
can verify, that IjlSfl with u(x) — y'x"o(a;), X > Oi solves 
((TT|l with a periodic potential U{x) — 2a{l — x)\uo{x)\^ if 



*uo(a;) solves (|T^ . Alternatively, from Eq. 



one can find a lattice potential U(x) for a given wave 
field it(x), subject to the requirement of regular behavior 
of Uxx/u: U{x) — — oj ~ 2au^. By any of the above 
methods, one can verify that the Jacobi elliptic func- 
tions pq(Kx, 777), where p=c,s or d, and q = n [more gen- 
erally one can consider 7ipq(a;) = (Apq'^(Ka;, 777.) + 5) 
with B > \A\, a and /3 being real integers], subject to 
proper choice of a, give the lattice potential 



U{x) — Uosn^{Kx,m), 



(41) 



where Uq is a constant and we returned to the parame- 
terization given by k and 777. (Although one can continue 
using Aj, the latter now loose the meaning they had in 
the context of the 1ST.) The respective solutions are sum- 
marized in the Table Hill 



TABLE III: Solutions of the type Upq(x,t) — \/Apq{Hix,m). 
The parameters in the second column have to be chosen to 
provide A > 0. 



pq 




A 








sn 


{K^m —- 


Uo/2)a 


-^2(1 + 


777) 


cn 


{Uo/2 - 




«:^(277i - 1) 




dn 


{Uo/2m 




n^{2 - m) - 


Uo/m 



If the nonlinearity depends on time, we again make use 
of the substitution leading to the equation [c.f. l|^] 



ivt + Vxx — Uasn^{Kx, m)v — 2a'|w|^7; = i^v. (42) 

As above the right hand side of Eq. (|42|1 is considered as 
a small perturbation: \^{t)\ <^ 1. Then, in accordance 
with the adiabatic approximation, the parameters k and 
777 will be considered as slow functions of time t, and 
respectively will be designated as k and 777. The equa- 
tions of their evolution are obtained in the same way as 
we obtained above the equations for the dependence of 
Ai,2 on time. More specifically, they are Eqs. (|24|1 . H26|l. 
where the period now coincides with the period of the 
periodic potential and L and N are expressed through k 
and in: L — L(ii, rh), N = N{k, in). The time-dependent 
parameters must satisfy initial conditions as follows 



ih{0) 



(43) 



where respective links between 777 and k are given in the 

Table nni 



Uxx ^ {lo + U{x) + 2(Tu'^)u, 



(40) 



B. BEC with a positive scattering length: sn-wave 



which allows one to construct a lattice potential U{x), 
in an explicit form, when a stationary periodic solution 



Let us start with adiabatic dynamics of a periodic 
matter wave in a BEC with a positive scattering length 
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(5- = 1) and with 2k^to > Uq, loaded in OL (|^ . Substi- 
tuting solution for the sn-wave (the first row of Tablc lITT|l : 



K2m- ^e-^'^'(i+")*sn(KX,m) 



(44) 



into (|42|) one obtains: 



L = 



— A^= ^ [K(m) - E(to)] . (45) 



Now one can distinguish two different situations, where 
maxima of the atomic density are located (i) in max- 
ima (Uq > 0) and (ii) in minima {Uq < 0) of the pe- 
riodic potential (examples are shown in FigQ (a) and 
(c), respectively). It is interesting to notice that, al- 
though the last situation admits the limit m — s- 0, 
where the distribution is transformed into a sine-wave: 
^'sn \/\Uo\/2e~^'^ *sin(Ka;), it does not allow transi- 
tion to the linear limit in the sense of a small amplitude, 
the latter being possible only if f/o > 0, corresponding to 
C/o/(2k2). 

For the next step we substitute (|45|l into (|24|) and (|26|l , 
what gives 

dk k{2k'^rh — Uq) 

X [E(m) - K(m)] [F,{m) - m'K{m)] , (46) 
dm 2mrh' {2k^Th — Uq 



-K(m) [E(m) - K{r 



Here standard notation m'(C) = 1 ^ ™(C) is used and 
A,n = t/o [E2(m) -m'K2(m)] 

-I- 2k^m^^[E{m)-K{rh)f -mK^{rh)Y 

The obtained system (I46II can be analyzed in the lim- 
iting cases. We start with the limit m — > at f/p < 0, by 
imposing the initial conditions 



to(0) = 0, 



2tt/L 



(47) 



(the case becomes trivial in absence of the periodic po- 
tential, when C/q = 0). Using asymptotics l|B2|) and con- 
dition TO ^ 1 we obtain from (|45|l 



= 1 



TO 



(48) 



Next, requiring 2k^to <C \Uq\^ what for k 1 is equiva- 
lent to |k— k| ^ Uq, we deduce from the first of equations 



k" 



1 



4|(7o|C 



16k2 



(49) 



where we used the condition 4J7oC <C 1, necessary for the 
smallness of to. 

In the opposite limit to ^ 1 one can use the asymp- 
totics (IBlll to obtain 



TO 



1 



-Lk/2 



(50) 



Now the potential amplitude does not affect the asymp- 
totics in the leading order. This is expectable behavior 
because it corresponds to large particle densities, when 
two body interactions dominate the effect of the periodic 
potential. 

In Fig^ we show examples of numerical solutions of 
equations (|46|l for three different situations with respect 
to Uq and intermediate values of the wave parameters. 
In all the cases the wave is transformed in a sequence of 
dark solitons. The external potential affects the velocity 
of the process. 




FIG. 1: Examples of sn-wave solutions at = (solid thick 
lines) and C, = C/in (solid thin lines) in the periodic potential 
ijlTjl (dashed lines) for (a) Uo = 0.05, and (fin = 5; (b) Uo = 
0, (fin = 6, and (c) Uo = —0.05, C/™ = 6. The inset shows 
dynamics of the parameters and ^2, which connection with 
m and k is given in Table HTl In all the cases m = 0.11 and 
K = 0.75. 



C. BEC with a negative scattering length: cn-wave 

Now we consider deformation of the periodic cn-wave 
resulting in creation of a train of bright matter solitons 
in a BEC with the negative scattering length, a = —1, 
loaded in OL H41(l . The cn-wave unaffected by FR is 
given by the second row of the Table ITTll Thus we look 
for a solution in the form 



i)-^«)*cn(Kx,TO) . (51) 



Notice that, as above, choice of the magnitude of the OL 
amplitude has the constrain Uq < 2k^'m. Unlike in the 
previous case, now the value of the amplitude of the OL 
changes not only the initial atomic distribution but also 
the frequency of the solution. 
Now we have 



L 



N : 



4K(to) 
k _ ' 

4k TO - 



2Uo 



Km 



[E(to) - to'K(to)] 



(52) 
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and equations of the adiabatic approximation (|24() . 1261) 
can be written down in the form: 



dk k{2k fh — Uq) 



[E(m) 



dm 2{2k^rh — Uo)rfirh' 
where 



K(m) [E(m) - m'K(m)] 



(53) 



+ 2k^m {E^(to) + m K(to) [K(m) - 2E(m)]} . 

As before we consider the hmiting cases starting with 
the initial conditions l|47|l what is possible only for neg- 
ative C/q- The link H48|) holds in the case at hand, while 
the dynamics of k is governed by the formula 



k" 



( 1 



4|^o|C 
16k2 - |[/n 



(54) 



The obtained result reveals an interesting feature: de- 
pending on the amplitude of the potential barrier be- 
tween two adjacent minima, Uq, the amplitude of the 
periodic wave is either increasing, if — IGk^ < J/q < 0, 
or decreasing, if Uq < —16k^. This is a general prop- 
erty, rather than the phenomenon observed in the limit 
of small m. Indeed, for a given Uq, the denominator Ad 
in (|53|l acquires zero value for some rh and k, while the 
numerators in H53|l are positive. In small vicinities of 
each of the curves, where Acn is close to zero the adia- 
batic approximation in the form H53|l is not valid. 

To explain the described behavior we notice that de- 
creasing of the wave amplitude is possible only for nega- 
tive Uf). In the case under consideration, where Uq < 0, 
atoms are collected around maxima of the potential 
(what happens due to attractive interactions and is im- 
possible in the linear limit). Thus the condensate under- 
goes the effect of two forces acting in opposite directions: 
one originated by attractive interaction among particles 
and another force, tempting to separate condensates, due 
to the lattice potential. When the height of the potential 
barrier becomes larger than some critical value the sec- 
ond force becomes dominant, and the condensate tends 
to "split" in the vicinity of each maxima resulting in de- 
crease of the amplitude and shift of the atomic density 
toward the potential minima (notice, that this situation 
reminds the case depicted in Fig^). 

Returning to the asymptotic (jSljl we observe that rh 
becomes negative if Uq < —16k^ and thus the atomic 
density is approximated by the function 



2 



k'^\iTi\e 



i(-K^(2tmH-l) + |C/o|)t 



cd vi + N 



1 + \m\ 



(55) 



When the amplitude of the periodic wave growth, i.e. 
when —16k^ < Uq < 2K^'m, one observes generation of 



the train of bright solitons. An example of numerical 
solutions of equations of the adiabatic approximation for 
different intensities of periodic potential is presented in 
FigH 



0.1 - 
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FIG. 2: Examples of cn-wave solutions at = (solid thick 
lines) and ( = (fin (soUd thin lines) in the periodic potential 
(dashed lines) for (a) Uo = 0.05, and (fin — 4.5; (b) Uo = 0, 
Qin = 4, and (c) Uo — —0.05, (fin = 3.5. The inset in (b) 
shows dynamics of the parameters 77 and ^ given by H53^ . In 
all cases m = 0.11 and k = 0.75. 



D. BEC with a negative scattering length, dn-wave 

As the last example of the adiabatic dynamics we con- 
sider evolution of the periodic dn-wave. As in the previ- 
ous case the attractive nature of the interaction (ct = —1) 
and application of the FR result in creation of a train of 
bright matter solitons. A solution of (|42|l is searched now 
in the form (see last row of Table IIII|I 




«(K^(2-m)-C/o/m)t^^(^^^^) (56) 



L = 



2K(m) 



N = 



2{2k'^m - C/o)E(m) 



K Km 
and equations of the adiabatic approximation read: 

dk k{2k^m — Uq)^ 

dC^ 



dm 2{2k^m — UQ)fkm' 



E(to) [E(m) - m'K(m)] 



rfC 
where 



Adn 



E(TO)K(m) 



(57) 



(58) 



Adn = C/o [E2(m) +m'K2(m)] 

+ 2k^rfi [E^(to) - mK^(m)] . 

An example of dynamics of the dn-wave, obtained by 
numerical solution of the system H58|) is shown in Fig|3| 
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FIG. 3: Examples of dn-wave solutions at ^ = (solid thick 
lines) and (" = Qin (solid thin lines) in the periodic potential 
(dashed lines) for (a) Uo = 0, and = 0.25; (b) Uo ~ 0.1, 
Cfin = 3. The inset in (a) shows the respective dynamics 
of parameters 771,2 given by l|33|l . In all cases m — 0.11 and 
K = 0.75. 

V. BEC IN AN OL SUBJECT TO FR. 
INCLUDING BOUNDARY EFFECTS 

In the last two sections we considered situations, where 
the parabolic trap is long enough in the axial direction 
and dynamics of the BEC can be treated within the 
framework of the adiabatic theory. In other two cases, 
where boundary effects are relevant [Eqs. ^ and (|13|l ] 
numerical analysis has to be employed. This rises the 
problem of a choice of adequate initial conditions (i.e. the 
initial conditions which in the absence of FR are station- 
ary solutions, such that the whole dynamics is originated 
by the FR only). Case 2, i.e. the NLS equation with 
an effective mass (O, was studied in In the present 
section we are concerned with the Case 4. 

To this end, we impose the following temporal behavior 
of the nonlinearity ^ 

g{t) = exp(VT) (59) 

where r is a constant. Then wave evolution is governed 
by the equation [c.f. ((T^ ] 

iil}t = -iixx + Uosn^{KX,m)ip + 2(7e*/^|'0pV' 

Wx^ip. (60) 

We also suppose that initial scattering length, as(0) 
is small enough, such that the two-body interactions 
are negligible. This implies smallness of the renormal- 
ized density, \tp\'^ <C 1. The aspect ratio is consid- 
ered to be ~ 10^'^ 10~^ what allows us to con- 
sider initial distributions in a form of modulated waves: 
'ip{x,t = 0) = exp{~h'x'^ /2)upq{x), close to the ground 
state in the absence of FR. 



A. Deformation of a cn-wave 

First we consider deformation of a periodic cn-wave re- 
sulting in creation of a train of bright solitons in a BEC 



with negative scattering length embedded in OL l|41(l . 
A typical example of numerical solution of Eq. H6U|1 with 
initial condition taken in a form of modulated cn-wave is 
presented in FigQl We show two subsequent regimes. 
First, on the time interval < t < 19 the FR was 
switched on and after that at 19 < t < 60 all poten- 
tials and FR were switched off, and the condensate was 
allowed to expand (freely in one dimension). 




FIG. 4; Dynamics of the cn-wave. The parameters are 1/ = 
0.0005, r = 2, K = 0.5, m = 0.1, and Uo = 0.05. In the inset 
growth of the number of solitons is shown. 

In Fig0| one can see that growth of the scattering 
length leads to sharpening of peaks of the distribution 
and to increase of their amplitudes. At an intermedi- 
ate moment of time (here t = 19) we obtain strongly 
pronounced train of pulses, which is naturally to asso- 
ciate with bright solitons. Strictly speaking the concept 
of a soliton is mathematically well defined in integrable 
systems, like NLS equation JT)), where each solution is 
associated with an eigenvalue of discrete spectrum of the 
ZS spectral problem (see Appendix A). In order to make 
use of this mathematical definition in our case, we refor- 
mulate it in the following way. 

Let us assume that a wave evolves during the time 
interval < t < tg. Then we assume that at to the trap 
and lattice potentials are switched off and the output 
pulse, i.e. one at t = to, is governed by the unperturbed 
NLS equation. Then we compute the number of NLS 
solitons generated by the wave profile given at i = to- To 
this end we substitute ^(a;,t = to) into the ZS spectral 
problem and study its spectrum. 

The result of this approach is shown in the inset of 
Fig0| At each time step {At = 0.005) in the interval 
< < < 19 the function ip{x,t) is considered as an ini- 
tial condition for the respective NLS equation. One can 
see that the initial cn-wave does not have solitons while 
increase of the nonlinear coefficient with time results in 
exponential growth of the number of solitons (to compute 
the number of solitons we used the numerical approach 
due to ,25J). 

In order to visualize the generated solitons in Fig^we 
show how the localized pulses, created due to the FR, 
propagate after all potentials are switched off at t = 19. 
In Fig^we can see that indeed the condensate consists of 
a set of solitary pulses propagating outwards the center 
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without significant changes of their shapes. This behav- 
ior supports the fact that we indeed generated train of 
bright sohtons. Meantime it should be mentioned that 
at i = 19 the number of sohtons calculated by solving 
spectral problem exceeds the number of solitons visible 
in Fig0] This can be explained by the facts that peaks in 
FigEl can present bound states consisting of several soli- 
tons which move with the same velocities and that small 
amplitude solitons are invisible on the scale of picture. 



B. Deformation of a dn-wave 

Now we consider deformation of the dn-wave resulting 
in creation of a train of bright matter solitons in a BEC 
with the negative scattering length. We solve Eq. H60() 
with initial condition taken as a dn-wave (third row in 
Table III) modulated by the Gaussian background. As 
in the previous case, increase of the strength of the inter- 
particle interactions by means of FR, results in growth 
and compression of the distribution picks, transforming 
them in a train of bright solitons (see the insert in FigEJ. 




FIG. 5: Dynamics of the dn-wave. The parameters are v = 
0.0005, T = 4, K = 1, m = 0.1, and (7o = 0.2. 

In contrast to the previous case, however, the gener- 
ated train of solitons is not robust, and is destroyed after 
some time of free propagation, corresponding in Fig[51to 
time interval t > to where the wave propagates freely. 
The difference in the behavior of trains of solitons can 
be understood if one takes into account that the dn-wave 
has a period of the lattice potential while cn-wave has a 
double period. As a result all bright solitons generated 
from the dn-wave are in-phase, while solitons generated 
from cn-wave have alternating phases (and thus the near- 
est neighbors are out-of-phase). 



C. Deformation of an sn-wave 

As the last example of this section we consider the gen- 
eration of dark solitons in a BECs with a positive scat- 
tering length and with initial condition taken in a form 
of the sn-wave (the first row of Table III), modulated 
by the Gaussian background. This case is very different 



from the two previous ones. The main reason is that 
now we are dealing with a periodic wave transforming in 
a train of dark solitons, which exist against a background. 
Either by change of the nonlinearity or by allowing free 
expansion of the condensate the background is strongly 
affected, what becomes an obstacle for numerical (and 
experimental) realization of the scenario described in the 
subsection IIII Dl It turns out however that it is possible 
to reduce the effect of the condensate expansion. 

To this end we first consider the NLS equation with 
constant nonlinear coefficient {g ~ I) 

idti^ = -dl^p + iy^x^TP + 2\TP\^ij, (61) 

and consider stationary solution tp « utf{x) e^*^* in the 
Thomas-Fermi approximation: |mtfP = (/^ — iy'^x'^)/2. 

Next we assume that FR is switched on and try to com- 
pensate growth of the nonlinearity by time dependence of 
the chemical potential /i(t). Then the respective solution 
is approximated hy ip — utf{x) e~'^(*) where M{t) is to 
be determined. In order to compensate the change of the 
nonlinear part we also vary in time the strength of the 
parabolic potential u — > v{t). Then the Thomas- Fermi 

approximation takes the form ImtfP = '^*^g(t)^^ • 
order to preserve the background profile during action of 
FR we require utf and utf to be equal what leads to 
the conclusion that the parameter ^^(i), defining change 
of the aspect ration of the cigar-shaped t rap, must change 
in time according to the law = v^J g(t). 

Thus, "stable" generation of a train of dark solitons 
can be achieved by simultaneous application of the FR 
and the respective trap potential what is described by 
the model 

+2e*/^|V'|V- (62) 

Examples of initial and final shapes of the condensate 
density affected by FR and adiabatically varying trap po- 
tential within the framework of the model l|62|l are shown 
in FigEl 



0.0001 




FIG. 6: Dynamics of dark solitons in the parabolic potential 
after switching off the FR and periodic potential at t = 21. 
Dashed, thin and thick lines correspond to t = 0; 21; 25. Pa- 
rameters are u = 0.0002, t = 2, k = 0.5, m — 0.1, and 
Uo = 0.05. 

To identify a train of dark solitons, avoiding free ex- 
pansion of the condensate, we switch off OL, FR and 
modulation of the trap potential, holding the last one 
unchanged. As one can see form the FigEl where the 
switching time is t = 21, the dark solitons start to move 
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toward the center of the condensate. This happens be- 
cause of reflection of soHtons by the potential boundaries 
[23 |. which due to introduced modulation of v the trap 
potential has become more narrow than it was at i = 0. 
The difference in the depth of each of solitons causes dif- 
ferent velocities and distance which solitons pass, what 
is also seen in FigEl 

VI. 3D NUMERICAL SIMULATIONS. THE 
GENERAL CASE 

As it is clear ID models considered above cannot ac- 
count all features of the BEC dynamics, what becomes 
especially relevant in the case of growing nonlincarity. In 
particular, the later can change significantly the back- 
ground resulting in failure of the low density approach. 
That is why we study numerically 3D GP equation 
for a cigar-shaped condensate having radially symmetric 
configuration. Using the same dimensionless variables as 
in Eq.Q we rewrite the GP equation in the form 

-I- C/oK2sn2(Kx,TO)*-|-87rcre*/^|«'p5' (63) 

where all tildes are suppressed and it is assumed that the 
nonlincarity varies according to (|59|) . 

By using the assumption about initially weak nonlin- 
carity one can consider initial profiles of the condensate 
as ID periodic solutions Upq(a:;), taken from the Table III, 
where must be substituted by k^U^ due to the nor- 
malization accepted in Sec. IIVI [see Ea. (|41|) ]. modulated 
by the 3D Gaussian envelope: 

|vI/(x,r^,<-0)|2 = ^i2^(x)e-'■i"''-^ (64) 

Such a distribution is close to the ground state and its 
temporal evolution occurs mainly due to change of the 
scattering length. 

Like in ID case we use implicit Crank-Nicolson scheme 
for numerical integration of the radially symmetric 3D 
GP equation. The scheme stability are fulfilled by con- 
trolling the ratio between time and space steps and the 
results are confirmed by using different time and space 
grids. On each time step of calculations in each site of the 
grid we check the convergence of the numerical scheme. 
In all 3D calculation we used grid of 400 x 400 points with 
the spatial step 0.2 and temporal step 0.001. 

A. Dynamics of cn- and dn-waves in a BEC with a 
negative scattering length 

First we consider adiabatic dynamics of a BEC with 
initial density distribution H64|l where Kpq(a;) is taken to 
be the cn-wave (see Table III). To be specific, the nu- 
merical results reported below in the physical units cor- 
respond to A/" « 5 X 10'^ atoms of ^Li loaded in an elon- 
gated trap with ojo = Stt x 1.2 Hz and = 27r x 57.5 Hz 



(the aspect ratio being v « 0.02). Taking into account 
that by changing the magnetic field near resonant point 
725 G one can change scattering length by order 10^-;- 10'^ 
the initial magnitude of the scattering length is cho- 
sen to be a<;(0) = —0.1 nm. This corresponds to transfer 
of the condensate from quasi-linear limit to the limit of 
relatively strong two-body interactions. The respective 
initial density distribution is shown in Fig [7^. 




FIG. 7: Dynamics of the condensate density, |»I'(r)|^, obtained 
by numerical solution of GP equation 16311 with attractive 
interactions [a = —1). Initial density is taken as 16411 where 
parameters of cn-wave are m = 0.1, k = 0.5 and v — 0.02. 
The amplitude of OL is Uo = 0.2Er with Er 0.7 peV. (a) 
Initial, (b) intermediate (to = 4.4 ms) and (c) final (tf = 
9 ms) density profiles. The characteristic time of FR r = 
1 ms. The size of each image is 180 x 40 /im. 

While FR is switched on (i.e. for t < 4.4 ms) the 
periodic density distribution displays sharpening in each 
potential well acquiring a shape of quasi-lD soliton train 
a,t t — 4.4 ms (Fig[7|D). As one can see from Figs[7^,b, 
appearing of localized excitations is accompanied by ap- 
preciable narrowing of the transverse profile of the con- 
densate (evidently not taken into account in the ID ap- 
proximation developed above). Change of the axial di- 
mension of the condensate is negligible. We also notice 
that due to choice of relatively weak initial nonlincarity 
collapse does not occur. 

Next, at the intermediate time to = 4.4 ms [what cor- 
responds to as{to) ~ 81.5as(0)] we switch off FR, OL 
and trap potential in the longitudinal direction (leaving 
only parabolic potential in the radial direction). This 
allows condensate to expand freely in the axial direc- 
tion (a kind of a waveguide configuration). A density 
distribution after 4.6 ms of free unidirectional expansion 
is shown in Fig[7|;. Now the transverse profile of the 
condensate changes weakly while localized pulses move 
outwards the center of the condensate. In the ID ap- 
proximation such pulses were identified as solitons. Now 
we also observe that low amplitude solitary pulses (they 
are located on the wings of the condensate) dispersion 
dominates nonlincarity and one observes spreading out 
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of each pulse. Large amplitude pulses preserve well pro- 
nounced spatially localized structure, while moving along 
the effectively created waveguide in the condensate. 




FIG. 8: Dynamics of the condensate density, |^'(r)p, obtained 
by numerical solution of GP equation 16311 with attractive 
interactions {a — —1). The initial density distribution l|64^ is 
taken as a dn-wave with m = 0.1, k = 1 and with f = 0.02. 
(a) Initial, (b) intermediate {to = 5.2 ms) and (c) final (tf — 
10 ms) density profiles of the condensate. The amplitude of 
OL is Uo = 0.2Er with Er 0.35 peV, and And r = 2 ms. 
The size of each image is 180 x 80 fim. 

Next we consider evolution of the dn-wave (where in 
the initial density distribution (|64|l Upq(x) = Udnix) from 
Table III) . Now we use the same magnetic trap as in the 
previous case, taking J\f w 10^ of ^Li atoms. Numerical 
simulations are carried out in the same way as this was 
done in the previous example. FR is switched off at time 
to ~ 5.2 ms (what corresponds to as(to) ~ 14as(0)) and 
free ID expansion of the condensate along the waveguide 
during 4.8ms between Iq and tf = 10 ms is considered. 

While the FR is switched on one observes strong 
sharpening of the distribution maxima (much more pro- 
nounced than in the case of the cn-wave shown in FiglTj) 
with significant narrowing of the transverse profile (see 
FiglHb)- Meantime, during the free expansion, the con- 
densate displays less regular structure of the train of soli- 
tary pulses (see FiglS]:), what was predicted by the ID 
theory and explained by the fact that now we have a 
sequence of the in-phase neighbor pulses. Another prop- 
erty, found numerically is that oscillatory behavior of the 



transverse distribution appears (the effect neglected by 
ID theory). 



B. Dynamics of an sn-wave in a BEC with a 
positive scattering length 

As a final example we consider the case of a BEC with 
repulsive interactions {a — I) where in the initial density 
distribution H64|l the solution Upq is taken in a form of the 
sn-wave (i.e. as Usn from Table III). In FiglHlwc present 
density profiles of the condensate of A/" ~ 5 x 10** atoms of 
*^Rb, confined by the trap with axial cjo = 27r x 0.11 Hz 
and radial uj± = 27rx 1.17 Hz harmonic oscillator frequen- 
cies, at different moments of time. It is supposed that 
initially two body interactions are negligible. The last 
condition for the given number of atoms can be achieved 
by proper choice of the external magnetic field. For ex- 
ample, according to 01 in the case of '^^Rb atoms the 
magnetic field about 170 G results in 0^(0) = 0.1 nm, 
what according to Q gives the small parameter e « 0.38, 
providing the required weakly nonlinear limit. Then by 
approaching the resonant point (w 156 G) one can change 
magnitude of scattering length by 10^ times. 

As before FR is switched on att — Q leading to the evo- 
lution shown for tg = 30 ms in Figs|HJi,b. One observes 
decreasing of the amplitude of the density and formation 
of a sequence of well pronounced excitations having the 
shape of a train of quasi-one-dimensional dark solitons. 

As before at some time, precisely at to = 30 ms 
when the scattering length achieves the value as(to) = 
148as(0), we switch off both FR and OL, but unlike in 
the previous cases, we leave both the transversal and the 
axial parabolic traps. The condensate continues to ex- 
pand in the transverse direction what leads to decrease 
its density and strong deviation from the original quasi- 
one-dimensional geometry. 



C. 3D dynamics with ID configuration. 
Comparison 

Numerical study performed in the previous subsections 
used elongated traps with typical for experiments values 
of the linear oscillator frequencies. While they illustrated 
qualitative agreement with prediction based on ID mod- 
els, they strictly speaking did not satisfy completely the 
conditions of the quasi- ID dynamics because of relatively 
wide transverse dimension. As a consequence 3D effects 
were not negligible especially at large nonlinearities. It 
is therefore a relevant problem to provide quantitative 
comparison of ID and 3D models. We will do this in de- 
tails for the 4-th Case listed in Table I, when the trans- 
verse size of the condensate is much smaller than the 
lattice constant. More specifically, we compare dynamics 
of the 3D condensate in the longitudinal direction gov- 
erned by (|63|l with dynamics of its ID counterpart de- 
scribed by Eq. H12|l . We use the dimensionless variables of 
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^^^^^^^ 

^^^^^^ 

(c) ^^^^^^^^^^ 

^^^^^^^ 



FIG. 9: Dynamics of the condensate density obtained by nu- 
merical solution of GP equation 1631 with repulsive interac- 
tions {a — 1). Initial density is taken as the sn-wave with 
m = 0.1 and k = 0.75 modulated by the Gaussian function 
with v — 0.1. (a) Initial, (b) intermediate {to = 30 ms) and 
(c) final {tf =45 ms) density profiles of the condensate. The 
OL amplitude is Uo « -0.2Er with Er « 0.034 peV, and 
r ~ 6 ms. The size of each image is 150 x 60 fim. 



3D problem, introduced in the last subsection, and take 
into account connections between ID and 3D variables, 
coming from the multiple-scale expansion: xijj ^ ex^u, 
tiD ^ e^tsD, *iD ^'an/e, ^ vs.d/^^- Natu- 
rally, for comparison we use the same initial distribu- 
tions in the longitudinal directions taken in a form of cn- 
[Figslini(a)-(c)] and dn- waves [Figs [101 (d)-(f)] for a BEC 
with attractive interactions and in a form of an sn-wave 
[Figs^l (g)-(i)] in the case of repulsive interactions. 



\A 










1 (e) 




(0 




FIG. 10: The cn-wave dynamics: (a) Initial, (b) intermedi- 
ate (to = 22) and (c) final {tf = 42) density profiles. The 
parameters are Uo ~ 0.03, m = 0.1, k = 0.75, v = 0.001, 
T — 5, and e ~ 0.071. The dn-wave dynamics: (d) Initial, 
(e) intermediate {to = 27.5) and (f) final {tf = 36) density 
profiles. The parameters are Uo = —0.1, m — 0.1, k — 0.75, 
v — 0.002, T = 10, and e ~ 0.068. The sn-wave dynamics: (g) 
Initial, (h) intermediate {to — 7) and (i) final {tf =9) density 
profiles. The parameters are Uo = —0.02, m — 0.001, k = 1, 
V = 0.002, r = 1, and e « 0.1. 

In FigllOlone observes a remarkable accuracy of the ID 
approximation, small deviations from which occur due to 



the transverse dynamics of the condensates, which be- 
comes more pronounced at large nonlinearities. 



VII. CONCLUSION 

To conclude we investigated formation of trains of 
bright and dark matter solitons in an elongated Bose- 
Einstein condensates embedded in optical lattices by 
means of adiabatic change of the value of the scattering 
length. For a strongly elongated cigar-shaped condensate 
we considered reductions of the original Gross-Pitaevskii 
equation to different one-dimensional models, depending 
on the parameters of the problem (the results are summa- 
rized in Table I). This was done in a controlled manner by 
using multiple-scale expansion method. Special attention 
has been paid to discussion of the link between periodic 
and solitary wave solutions of the one-dimensional non- 
linear Schrodinger equation, what justified an idea how 
to affect a periodic matter wave loaded in a lattice, in 
order to create a train of solitons. In a number of cases, 
where boundary effects are negligible, the phenomenon 
can be described analytically using the adiabatic approxi- 
mation. When both parabolic and periodic potentials are 
included explicitly in the one-dimensional model we had 
to involve numerical simulations in order to describe wave 
evolution. The consideration was concentrated on the 
three types of waves for which analytical expressions are 
available: cn-, dn- and sn-waves. Finally, the process of 
creation of the trains of bright and dark solitons by means 
of varying scattering length using Feshbach resonance 
in the presence of parabolic and periodic potentials was 
studied numerically withing the framework of the radi- 
ally symmetric three-dimensional Gross-Pitaevskii equa- 
tion. It has been shown that for the most typical experi- 
mental parameters one-dimensional approximation gives 
relatively good qualitative description of the dynamics. 
Meantime, it turned out that full three-dimensional con- 
sideration may be necessary because of significant de- 
formation of the transverse structure of the condensate, 
especially in the case of dark solitons, which are very sen- 
sitive to the amplitude and structure of the background, 
affected by the Feshbach resonance. 

If all conditions of the applicability of the multiple scale 
expansion are satisfied, one observe remarkable quantita- 
tive agreement between the real three-dimensional mod- 
els and it one-dimensional counterpart. 
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APPENDIX A: SOME ELEMENTS OF THE 
INVERSE SCATTERING TECHNIQUE 

The starting point of the 1ST is the fact that Eq. l(TI| 
appears as a compatibility condition for two systems of 
differential equations (see e.g. pol |') 



l/>. = U|/), |/)t-V|/) 
where |/) = col(/i,/2), fi.2 = fi.2{^,t), 

U 



( 2i 



V = 



2i 



(Al) 

(A2) 
(A3) 



and A is called a spectral parameter. The first of systems 
(|A1I) is referred to as the Zakharov-Shabat (ZS) spectral 
problem. 

Let us introduce (/| = (— /2,/i) and define the inner 
product {f\g) = /152 — /251 and a projector matrix 



l/>(5l 



-/1.92 

-/2.g2 /251 



(A4) 



If I/) and \g) are two independent solutions of each 
of the systems (|A1|) . then {f\g) is a Wronskian of each 
of them, and hence it does not depend neither on x nor 
on i.e. {f\g) = p(A), where p(A) is a function of the 
spectral parameter. Let |/) = col(/i, /2) be a solution of 
(|A1I) with a given ?/;. Then the symmetry properties of 
the matrices U and V imply that \g) — col(/2,(7/i) is a 
solution of (|A1|) . as well. Using the pair of the solutions 
I/) and \g) one obtains that {f\g) = 5-|/ip-|/2p =p(A). 

Next, one ensures that P solves equations 



= [U,P], P* = [V,P], 
from which it follows that 

AP2 _ 4P12P21 - P{\) 



(A5) 



(A6) 



where AP — P22 — Pii, Pij being entries of the matrix 
P, and P(A) is a function of A only. For the chosen 
pair of solutions one verifies that P21 = ^Pi2- Thus 
^12^21 = ^|/i|'|/2p and AP ^ + I/2I2. Conse- 

quently, P(A) = p^(A) is a real positive function of A. 

Noticing that U and V are matrix polynomials in A, 
one can conjecture that the simplest nontrivial solution 
I/), and hence P, should also be searched in forms of 
polynomials with respect to A. As far as V is a quadratic 
(matrix) polynomial of A, one has to consider |/) to be 
a second (or higher) order polynomial with respect to A. 
Since P is quadratic with respect to the elements of |/) 
one concludes that the simplest nontrivial form of P(A) 
must also be a polynomial of the forth degree of A: 



P(A) 



(A 
A^ 



Ai)(A 



-A2)(A-A3)(A 
C2A^ + ciA 



A4) 



Co- 



(A7) 



Here Xj{j = 1,...,4) are given complex numbers, and 
Cj are trivially expressed through Xj. Then, = ReAj 
and rjj = ImXj can be considered as parameters char- 
acterizing respective solutions of the nonlinear problem 
(|14|1 . The requirement for P(A) to be real imposes four 
constrains on the parameters Xj. Indeed, considering 
P(A) one concludes that either (i) all Xj are real, i.e. 
rjj — and S^j parametrize the problem, or (ii) two of 
them are real and two ones are complex conjugated, say 
ImAi=ImA2 = and A3 = A4, then the parameters are 
{Cij ^2, Csi or (iii) there are two pairs of cojnplex con- 
jugated parameters, say Ai = A2 and A3 — A4, i.e. the 
parameters are {Ci, 61 J?!; 

Formulas (|A6|I and \K1\ suggest that P^ should also 
be searched in a form of polynomials with respect to the 
parameter A: 



AP = A^ + AAi + Ao, 



P11 



AP 



(1) 
12 



P 



(0) 



12 



where explicit form of the coefficients Aj and P^-i are to 



(A8) 



be found. Comparing powers of A in (|A6|) we obtain 



Al = 



C3 



2 ' 

C3A0 - 4a (Pi^2 



Ao = 2a|Pi2 



(1)|2 



1 



C2 



P 



(0) 



p(0) p(l) 
^12 ^12 



Cl, 



A^-4a|P| 



(0)|2 



Co- 



,(A9) 

(AlO) 
(All) 



From the first of Eqs. ljA5|l we compute the equation for 
9xPi2 whose consistency with representations (|A8|) and 
(|A9(I results in the following connecting formulas 



\J a 



1 



(0) 



i'. = ^Pk 



2^3, 



(A12) 

(A13) 
(A14) 



Comparing IIA14p with {Tj) one verifies that the two equa- 
tions coincide subject to the conditions C2 — 2lu and 
C3 = what is the same as 







Al + A2 + A3 + A4 



(A15) 



(A16) 



It follows from (|A12p that P^2^ is a solution of ifTHl . 
Taking into account the particular form of solution l|15() 
one immediately concludes from (|A10|) . I|A12(I . I)A13(I . 
and (|A16p . that ci = 0. Finally, comparing (|A11|) 
with 1161) and taking into account H15() one obtaines 
C = jg{c2 — 4co) what is the same as 



C= -(C02-A1A2A3A4) 



(A17) 



The last condition leaves us with only two independent 
parameters. 
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A convenient representation of the results is given 
by the locations of the parameters Xj on the Riemann 
surface of the function p(A), which is shown in Fig llll 
(Figsn](a), (c), and (e) correspond to the above cases 
(i), (ii), and (iii)). 



(a) 



(o) 



(b) 



:i 








4. f 




t J 


! -in. 


(b) 




- € + in 




X — 


X 






-s-in 


X 



>:177 



y.-V\ 



-177 



If? 



■XT) 



FIG. 11: The Riemann surface of the function p(A) = \J -P(A) 
consists of two complex plains, one of which is shown in Figs, 
(a) for sn-wave: Ai = —A3 = ^1, A2 — — A4 — ^2; (c) for 
dn-wave: Ai = A3 = ir\\, A2 = A4 = 17^2; and (e) for cn- 
wave: Ai = —A3 = — A2 = — A4 = ^ + 177. The right panels 
show the modification of the Riemann surface subject to the 
limiting transition m 1 shown in the last column of the 
Table ^ The lines connecting the parameters Aj show cuts 
of the respective two-sheet Reamann surfaces. 



the ZS spectral problem, which determines a static dark 
soliton, while points ±^1 = ±2/3 are transformed into 
edges of the continuum spectrum. 

In the case of the NLS equation with a = —\ and zero 
boundary conditions, lima;_,-|-oo IV'I = 0, bright soliton 
solutions are determined by discrete spectrum located in 
the upper half-plane of the spectral parameter. If a soli- 
ton is static, the only discrete eigenvalue is placed on the 
imaginary axis, and is designated as ir\ fFigslllb and d). 
One can arrive at such situation by shifting the branch- 
ing points of corresponding dn-wave as it is shown on the 
passage between FigstTTb and d. The points 771,2 collapse 
into a point r\ corresponding to the discrete eigenvalue of 
the ZS spectral problem. Alternatively, one can arrive at 
a static bright soliton solution, starting with a cn-wave, 
and providing displacements of the branching points as 
it is shown by the arrow between Figslllb and f. 



APPENDIX B: SOME FORMULAS FROM 
ELLIPTIC INTEGRALS 

Throughout the text we use the asymptotic of the el- 
liptic integrals (see e. g. [l7Ll2^ ') when m — > 1: 



K(7ri) = In ■ 



[l + 0(l-m)]. 



E(7n) = 1 + O [(1 - 77i) ln(l - m) 



(Bl) 



and when m ^ 



K(m) = - 



E(77T.) 



1 + J + ( ^ ) + 0(to3^ 



3 \ m 



(B2) 



In order to prove that ^2 described by (|28|l is a mono- 
tonically decreasing function, we observe that K(m) > 
E{m), and show that 



Let us now assume that by some means one can pro- 
vide shifts of the spectral parameters in the complex 
plane. More specifically, we consider three situations 
corresponding to collapse of two pairs of the branching 
points, Xj, as it is shown by the arrows in Figllll 

At CT — 1, subject to transition ^2 — *■ and ^1 — > 
^ — 2p fFiglTTb) the Riemann surface of p{X) is trans- 
formed in the Riemann surface of the spectral parameter 
of the ZS problem, subject to the boundary conditions 
\\mx^±oo 4' = pe^^^ , where p and ip are real constants. 
Then the passage between (a) and (b) illustrates defor- 
mation of the branching points ^1.2 resulting in transfor- 
mation of the sn-wave into a single dark soliton. The 
point ^2 = corresponding to the discrete spectrum of 



2^2K{m) - (a + ^2)E{m) < 0, (B3) 

where m = (^1 — ^2)^/(^1 + ^2)^- To this end we 
introduce the elliptic modulus k — ^/m and rewrite 
(|B3|) as J^i(fc) = (1 - k)K{P) - E{k^) < 0. In or- 
der to prove the last formula we use that its left hand 
side is zero at A; = 0, and its derivative is given by 
dJ'i{k)ldk = -[E(fc2) + K(fc2)]/(1 + /c) < (for the 
last formulas we used j^^: dE/dk — (E — K)/fc and 
dK/dk = [E(fc2) - A:'2K(fc2)]/(fcfc'2) where k'^ + P = 1. 

By analogy, in order to show that 772 in (|33|) is mono- 
tonically increasing function it is enough to verify that 
dT2{k)/dk = -[K(fc2) _ E(fc2)][i _ k']/{kk') < where 
.F2(fc) = A:'K(P)-E(P)<0. 
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